// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                            avtImageFileFormat.C                           //
// ************************************************************************* //

#include <avtImageFileFormat.h>

#include <string>
#include <vector>
#include <visitstream.h>
#include <visit-config.h>

#include <vtkCellData.h>
#include <vtkCellType.h>
#include <vtkFloatArray.h>
#include <vtkPointData.h>
#include <vtkRectilinearGrid.h>
#include <vtkUnsignedCharArray.h>
#include <vtkUnstructuredGrid.h>

#ifndef MDSERVER
#include <vtkExtractVOI.h>
#include <vtkPNMReader.h>
#include <vtkPNGReader.h>
#include <vtkJPEGReader.h>
#include <vtkTIFFReader.h>
#include <vtkBMPReader.h>
#include <vtkStimulateReader.h>
#endif

#include <avtCallback.h>
#include <avtDatabaseMetaData.h>
#include <avtDataSelection.h>
#include <avtGhostData.h>
#include <avtLogicalSelection.h>
#include <avtParallel.h>
#include <avtSpatialBoxSelection.h>

#include <DebugStream.h>
#include <BadIndexException.h>
#include <InvalidFilesException.h>
#include <InvalidVariableException.h>

#define EXPOSE_DEPTH
#ifdef EXPOSE_DEPTH
#include <vtkZLibDataCompressor.h>
#include <FileFunctions.h>
#endif
using     std::vector;
using     std::string;

// ****************************************************************************
//  Method: avtImage constructor
//
//  Programmer: Chris Wojtan
//  Creation:   Thu Jun 3 09:50:31 PDT 2004
//
//  Modifications: 
//
//     Chris Wojtan, Mon Jun 21 15:36 PDT 2004
//     Moved most of the Image loading code into ReadInImage() instead of the
//     constructor
//
//     Mark C. Miller, Tue Nov  9 13:41:33 PST 2004
//     Moved code to load variable names from PopulateDatabaseMetaData to here
//     Moved code to determine file extension from ReadInImage to here
//
//     Hank Childs, Fri Mar 18 10:11:36 PST 2005
//     Move heavy lifting to Initialize method.
//
// ****************************************************************************

avtImageFileFormat::avtImageFileFormat(const char *filename, const DBOptionsAttributes *)
    : avtSTSDFileFormat(filename)
{
    fname = filename;
    image = NULL;
    haveReadWholeImage = false;
    haveImageVolume = false;
    haveInitialized = false;
    indexOfImageAlreadyRead = -1;
    indexOfImageToRead = -1;
}

// ***************************************************************************
//  Method: avtImageFileFormat destructor
//
//  Programmer: Chris Wojtan
//  Creation:   June 6, 2004
//
//  Modifications:
//
//    Mark C. Miller, Tue Nov  9 13:41:33 PST 2004
//    Replaced duplicate code to free things up with a call to
//    FreeUpResources
//
// **************************************************************************

avtImageFileFormat::~avtImageFileFormat()
{
    FreeUpResources();
    cellvarnames.clear();
}


// ****************************************************************************
//  Method: avtImageFileFormat::Initialize
//
//  Purpose:
//      Initializes the reader.  This is not done in the constructor, because
//      many of these objects may be instantiated at one time.
//
//  Programmer: Hank Childs
//  Creation:   March 18, 2005
//
//  Modifications:
//
//     Hank Childs, Thu Mar 17 15:33:44 PST 2005
//     If we are using the stimulate image format, don't add color channels.
//     Also add support for image volumes, where we pull out the extension.
//
//     Jeremy Meredith, Thu Aug  7 16:12:17 EDT 2008
//     Use const char* for string literals.
//
//     Jeremy Meredith, Thu Jan  7 12:52:02 EST 2010
//     Make sure dot position exists before using it; otherwise we crash.
//
//     Jeremy Meredith, Thu Jan 14 12:04:27 EST 2010
//     Oops, previous fix was neither correct nor sufficient.
//     Used modern for loop scoping to clarify.
//
// ****************************************************************************

void
avtImageFileFormat::Initialize(void)
{
    if (haveInitialized)
        return;
    haveInitialized = true;

    // find the file extension
    size_t start = fname.size();
    for(size_t i=0; i<fname.size(); i++)
        if(fname[i] == '.')
            start = i;
    if (start>=fname.size())
        fext = "";
    else
        fext = string(fname, start+1, fname.size()-1);

    debug4 << "File extension: " << fext.c_str() << endl;

    // If we have an image volume, do some special processing.
    if (fext == "imgvol")
    {
        ReadImageVolumeHeader();
        for(size_t i=0; i<subImages[0].size(); i++)
            if(subImages[0][i] == '.')
                start = i;
        image_fext = string(subImages[0], start+1, subImages[0].size()-1);
    }
    else
    {
        image_fext = fext;
        subImages.push_back(fname);
    }

    // Check to see if we have color channels.
    bool hasColorChannels = true;
    int numExtsWithoutColorChannels = 4;
    const char *extsWithoutColorChannels[4] = { "spr", "SPR", "sdt", "SDT" };
    for (int i = 0 ; i < numExtsWithoutColorChannels ; i++)
        if (image_fext == extsWithoutColorChannels[i])
            hasColorChannels = false;

    // load variable names
    if (hasColorChannels)
    {
        cellvarnames.push_back("red");
        cellvarnames.push_back("green");
        cellvarnames.push_back("blue");
        cellvarnames.push_back("alpha");
    }
    cellvarnames.push_back("intensity");
}


// ****************************************************************************
//  Method: avtImageFileFormat::ActivateTimestep
//
//  Purpose:
//      Activates this file format, allowing some initialization work to
//      happen.
//
//  Programmer: Hank Childs
//  Creation:   March 18, 2005
//
// ****************************************************************************

void
avtImageFileFormat::ActivateTimestep(void)
{
    Initialize();
}


// ****************************************************************************
//  Method: avtImageFileFormat::ReadImageVolumeHeader
//
//  Purpose:
//      Reads in a header file that lists a series of images that should be
//      combined as a volume.
//
//  Programmer: Hank Childs
//  Creation:   March 18, 2005
//
//  Modifications:
//
//    Hank Childs, Fri May 20 10:33:16 PDT 2005
//    Add support for relative paths.
//
//    Hank Childs, Mon Jun  6 16:19:23 PDT 2005
//    Read in the header on processor 0 only and broadcast to the other
//    processors.
//
//    Brad Whitlock, Thu May 11 16:22:03 PST 2006
//    Added code to strip off extra windows end of line characters from the
//    end of the file extensions so the reader is not confused later.
//
//    Mark C. Miller, Fri Apr 23 17:40:30 PDT 2010
//    Added warning about wasting processors.
// ****************************************************************************

void
avtImageFileFormat::ReadImageVolumeHeader(void)
{
    int rank = PAR_Rank();
    int size = PAR_Size();
    if (PAR_UIProcess())
    {
        ifstream ifile(fname.c_str());
        if (ifile.fail())
        {
            debug1 << "Cannot open image volume file \"" << fname.c_str() 
                   << "\"." << endl;
            int success = 0;
            BroadcastInt(success);
            EXCEPTION1(InvalidFilesException, fname.c_str());
        }

        specifiedZStart = 0;
        zStart = 0;
        specifiedZStep = 0;
        zStep = 1;
    
        const char   *filename = fname.c_str();
        char          dir[1024];
        const char   *p = filename, *q = NULL;
        while ((q = strstr(p, VISIT_SLASH_STRING)) != NULL)
        {
            p = q+1;
        }
        strncpy(dir, filename, p-filename);
        dir[p-filename] = '\0';
    
        char line[1024];
        while (!ifile.eof())
        {
            ifile.getline(line, 1024);
            if (line[0] == '\0')
                continue;
            if (strncmp(line, "Z_START:", strlen("Z_START:")) == 0)
            {
                zStart = atof(line + strlen("Z_START:"));
                specifiedZStart = true;
            }    
            else if (strncmp(line, "Z_STEP:", strlen("Z_STEP:")) == 0)
            {
                zStep = atof(line + strlen("Z_STEP:"));
                specifiedZStep = true;
            }    
            else
            {
                char line_with_dir[1024];
                if (line[0] == VISIT_SLASH_CHAR)
                {
                    strcpy(line_with_dir, line);
                }
                else
                {
                    sprintf(line_with_dir, "%s%s", dir, line);
                }

                // Trim any weird control characters off of the end.
                int len = (int)strlen(line_with_dir);
                char *cptr = line_with_dir + len - 1;
                while(cptr >= line_with_dir)
                {
                    if(*cptr < ' ')
                        *cptr-- = '\0';
                    else
                        break;
                }

                subImages.push_back(line_with_dir);
                debug4 << "Adding \"" << line_with_dir
                       << "\" to the list of subImages." << endl;
            }
        }

        //
        // Now broadcast what we've read to other processors.
        //
        int success = 1;
        BroadcastInt(success);
        int iTmp = specifiedZStart;
        BroadcastInt(iTmp);
        double dTmp = zStart;
        BroadcastDouble(dTmp);
        iTmp = specifiedZStep;
        BroadcastInt(iTmp);
        dTmp = zStep;
        BroadcastDouble(dTmp);
        BroadcastStringVector(subImages, rank);

        //
        // Issue warning of possible mistake in using too many processors.
        // Maximum parallelism is one MPI task per image.
        //
        static bool haveIssuedProcWarning = false;
        if ((size_t)size > subImages.size() && !haveIssuedProcWarning)
        {
            char msg[1024];
            snprintf(msg, sizeof(msg),
                "Because your imgvol dataset named \"%s\",\n"
                "has %d slices, it can be decomposed for parallel across at most %d processors.\n"
                "You are using %d processors and so %d of these are not adding any additional\n"
                "parallel speedup to operations involving it.",
                fname.c_str(), (int)subImages.size(), (int)subImages.size(), size, (int)(size-subImages.size())); 
            if (!avtCallback::IssueWarning(msg))
                cerr << msg << endl;
            haveIssuedProcWarning = true;
        }
    }
    else
    {
        int success = 0;
        BroadcastInt(success);
        if (!success)
        {
            debug1 << "Processor 0 could not open file \"" << fname.c_str() 
                   << "\"." << endl;
            EXCEPTION1(InvalidFilesException, fname.c_str());
        }
        int iTmp;
        BroadcastInt(iTmp);
        specifiedZStart = (bool) iTmp;
        double dTmp;
        BroadcastDouble(dTmp);
        zStart = dTmp;
        BroadcastInt(iTmp);
        specifiedZStep = (bool) iTmp;
        BroadcastDouble(dTmp);
        zStep = dTmp;
        BroadcastStringVector(subImages, rank);
    }

    haveImageVolume = true;
}


// ***************************************************************************
//  Method: avtImageFileFormat::CreateCacheNameIncludingSelections 
//
//  Purpose:
//      The image reader will return different data sets based on data
//      selections.  This method gives a description of what is returned
//      so that the resulting data can be reliably cached.
//
//  Programmer: Hank Childs
//  Creation:   December 20, 2011
//
// **************************************************************************

std::string
avtImageFileFormat::CreateCacheNameIncludingSelections(std::string s)
{
    int xmin, xmax, ymin, ymax;
    bool haveSelections = ProcessDataSelections(&xmin, &xmax, &ymin, &ymax);
    if (!haveSelections)
        return s;

    char str[1024];
    strcpy(str, s.c_str());
    size_t amt = strlen(str);
    for (size_t i = 0 ; i < selList.size() ; i++)
    {
        if ((*selsApplied)[i])
        {
            std::string s = selList[i]->DescriptionString();
            snprintf(str+amt, 1024-amt, "_%s", s.c_str());
            amt += strlen(str);
        }
    }
    return std::string(str);
}


// ***************************************************************************
//  Method: RegisterDataSelections 
//
//  Programmer: Mark C. Miller 
//  Creation:   November 3, 2004 
//
// **************************************************************************

void
avtImageFileFormat::RegisterDataSelections(
    const vector<avtDataSelection_p> &sels,
    vector<bool> *selectionsApplied)
{
    selList     = sels;
    selsApplied = selectionsApplied;
}

// ***************************************************************************
//  Method: ProcessDataSelections 
//
//  Purpose: Walk through all the data selections, decide which we can
//  service here and compose all that we can service into a single, logical
//  selection returned in the argument.
//
//  Programmer: Mark C. Miller 
//  Creation:   November 3, 2004 
//
//  Modifications:
//
//    Mark C. Miller, Tue Nov  9 13:41:33 PST 2004
//    Made it not process any selections for certain image file formats
//
//    Hank Childs, Thu Mar 17 15:33:44 PST 2005
//    Don't process data selections for Stimulate images.
//
//    Hank Childs, Tue Dec 20 16:52:07 PST 2011
//    Fix bug with logical selection check of 2D images.
//
// **************************************************************************

bool
avtImageFileFormat::ProcessDataSelections(int *xmin, int *xmax,
                                          int *ymin, int *ymax)
{
    bool retval = false;

    // some image file formats don't support selection on read.
    if ((image_fext == "png")||(image_fext == "PNG")||(image_fext == "tif")||
        (image_fext == "tiff")||(image_fext == "TIF")||(image_fext == "TIFF")||
        (image_fext == "jpg")||(image_fext == "jpeg")||(image_fext == "JPG")||
        (image_fext == "JPEG") || (image_fext == "spr") || (image_fext == "SPR")
        || (image_fext == "sdt") || (image_fext == "SDT") 
        || (image_fext == "imgvol"))
    {
        for (size_t i = 0; i < selList.size(); i++)
            (*selsApplied)[i] = false;
        return retval;
    }

    avtLogicalSelection composedSel;
    for (size_t i = 0; i < selList.size(); i++)
    {
        if (string(selList[i]->GetType()) == "Logical Data Selection")
        {
            avtLogicalSelection *sel = (avtLogicalSelection *) *(selList[i]);
            int strides[3];
            sel->GetStrides(strides);

            // currently handles only stride 1 
            bool stridesOk = true;
            int numDims = 0;
            sel->GetNDims(numDims);
            for (int j = 0; j < numDims; j++)
            {
                if (strides[j] != 1)
                {
                    stridesOk = false;
                    break;
                }
            }

            if (stridesOk == false)
            {
                (*selsApplied)[i] = false;
                continue;
            }

            // overrwrite method-scope arrays with the new indexing
            composedSel.Compose(*sel);
            (*selsApplied)[i] = true;
            retval = true;

        }
        else if (string(selList[i]->GetType()) == "Spatial Box Data Selection")
        {
            avtSpatialBoxSelection *sel = (avtSpatialBoxSelection *) *(selList[i]);

            double mins[3], maxs[3];
            sel->GetMins(mins);
            sel->GetMaxs(maxs);
            avtSpatialBoxSelection::InclusionMode imode =
                sel->GetInclusionMode();

            // we won't handle clipping of zones here
            if ((imode != avtSpatialBoxSelection::Whole) &&
                (imode != avtSpatialBoxSelection::Partial))
            {
                (*selsApplied)[i] = false;
                continue;
            }

            int imins[3], imaxs[3];
            for (int j = 0; j < 3; j++)
            {
                int imin = (int) mins[j];
                if (((double) imin < mins[j]) &&
                    (imode == avtSpatialBoxSelection::Whole))
                    imin++;
                
                int imax = (int) maxs[j];
                if (((double) imax < maxs[j]) &&
                    (imode == avtSpatialBoxSelection::Partial))
                    imax++;

                imins[j] = imin;
                imaxs[j] = imax;
            }

            avtLogicalSelection newSel;
            newSel.SetStarts(imins);
            newSel.SetStops(imaxs);

            composedSel.Compose(newSel);
            (*selsApplied)[i] = true;
            retval = true;
        }
        else
        {
            // indicate we won't handle this selection
            (*selsApplied)[i] = false;
        }
    }

    int mins[3], maxs[3];
    composedSel.GetStarts(mins);
    composedSel.GetStops(maxs);

    *xmin = mins[0]; 
    *xmax = maxs[0]; 
    *ymin = mins[1]; 
    *ymax = maxs[1]; 

    return retval;
}

// *****************************************************************************
//  Method: avtImageVileFormat::ReadInImage
//
//  Purpose:
//      Reads in the image.
//
//  Programmer: Chris Wojtan
//  Creation:   June 6, 2004
//
//  Modifications:
//
//    Chris Wojtan Jun 21, 15:38 PDT 2004
//    Moved most of the loading code into this function in order to speed up
//    the constructor
//
//    Jeremy Meredith, Wed Jul  7 13:49:18 PDT 2004
//    The MDServer doesn't need to read the images; it is making assumptions
//    about what is in the file already for speed purposes, and this lets us
//    get away with not using the extra VTK libraries.  Also commented out
//    the #include files at the top of the file.  Also, there was no case
//    for monochrome images, so I added one, and I threw an exception on an
//    unexpected number of channels.  And, made intensity be the numerical
//    average of the RGB channels instead of the sum.
//
//    Mark C. Miller, Thu Nov  4 18:21:25 PST 2004
//    Add some support for data selections. Not all formats can perform
//    selection during read. Those that can do. Those that cannot read the
//    whole image and then extract the region of interest after the read.
//    This is still useful to VisIt as it can reduce the amount of data
//    the engine keeps around
//
//    Mark C. Miller, Tue Nov  9 13:41:33 PST 2004
//    Backed off on the above change, a bit. Now, if a format doesn't support
//    selection on read, then it will read and keep around the whole image
//    If selection on read is supported but the selection is such that it winds
//    up reading the whole image, then it will read and keep the whole image.
//    Otherwise, it will read and keep only the portion requested.
//    Also, removed all code to allocate large, float vectors and populate them
//    to the GetVar call which is when they are actually needed.
//
//    Hank Childs, Fri Mar 11 10:05:51 PST 2005
//    Fix memory leak.
//
//    Hank Childs, Thu Mar 17 15:33:44 PST 2005
//    Add support for the stimulate image format.
//
//    Kathleen Bonnell, Wed May 17 14:14:04 PDT 2006 
//    Remove call to SetSource(NULL), as it now removes information necessary
//    for the dataset. 
//
//    Eric Brugger, Mon Nov 23 12:53:21 PST 2020
//    Correct the reading of Stimulate files to handle the case of a negative
//    step. Previously, a negative step would have resulted in monotonically
//    decreasing coordindates, which is invalid for a rectilinear grid.
//
// *****************************************************************************

void avtImageFileFormat::ReadInImage(void)
{
#ifdef MDSERVER
    return;
#else

    // if we've already read the entire image, then do nothing
    if (haveReadWholeImage &&
        ((indexOfImageToRead < 0) ||
         (indexOfImageToRead == indexOfImageAlreadyRead)))
    {
        return;
    }

    if (image != NULL)
    {
        image->Delete();
        image = NULL;
    }

    // process any data selections we can handle, here
    int xmin, xmax, ymin, ymax;
    bool haveSelections = ProcessDataSelections(&xmin, &xmax, &ymin, &ymax);

    haveReadWholeImage = true;
    indexOfImageAlreadyRead = indexOfImageToRead;
    if (haveSelections)
        haveReadWholeImage = false;

    xStart = yStart = 0;
    xStep = yStep = 1;

    int idx = indexOfImageToRead;
    idx = (idx < 0 ? 0 : idx);

    debug4 << "avtImageFileFormat::ReadInImage: image_fext=" << image_fext << endl;

    // select the appropriate reader for the file extension
    if ((image_fext == "pnm") || (image_fext == "PNM") || 
        (image_fext == "ppm") || (image_fext == "PPM"))
    {
        vtkPNMReader *reader = vtkPNMReader::New();
        if (haveSelections)
            reader->SetDataVOI(xmin,xmax,ymin,ymax,0,0);
        reader->SetFileName(subImages[idx].c_str());
        reader->Update();
        image = reader->GetOutput();
        image->Register(NULL);
        //image->SetSource(NULL);
        reader->Delete();
    }
    else if ((image_fext == "png") || (image_fext == "PNG"))
    {
        vtkPNGReader *reader = vtkPNGReader::New();
        reader->SetFileName(subImages[idx].c_str());
        reader->Update();
        image = reader->GetOutput();
        image->Register(NULL);
        //image->SetSource(NULL);
        reader->Delete();
    }
    else if ((image_fext == "jpg") || (image_fext == "jpeg") || 
             (image_fext == "JPG") || (image_fext == "JPEG"))
    {
        vtkJPEGReader *reader = vtkJPEGReader::New();
        reader->SetFileName(subImages[idx].c_str());
        reader->Update();
        image = reader->GetOutput();
        image->Register(NULL);
        //image->SetSource(NULL);
        reader->Delete();
    }
    else if ((image_fext == "tif") || (image_fext == "tiff") || 
             (image_fext == "TIF") || (image_fext == "TIFF"))
    {
        vtkTIFFReader *reader = vtkTIFFReader::New();
        reader->SetFileName(subImages[idx].c_str());
        reader->Update();
        image = reader->GetOutput();
        image->Register(NULL);
        //image->SetSource(NULL);
        reader->Delete();
    }
    else if ((image_fext == "bmp") || (image_fext == "BMP"))
    {
        vtkBMPReader *reader = vtkBMPReader::New();
        if (haveSelections)
            reader->SetDataVOI(xmin,xmax,ymin,ymax,0,0);
        reader->SetFileName(subImages[idx].c_str());
        reader->Update();
        image = reader->GetOutput();
        image->Register(NULL);
        //image->SetSource(NULL);
        reader->Delete();
    }
    else if ((image_fext == "spr") || (image_fext == "SPR") || 
             (image_fext == "sdt") || (image_fext == "SDT"))
    {
        vtkStimulateReader *reader = vtkStimulateReader::New();
        reader->SetFileName(subImages[idx].c_str());
        reader->Update();
        image = reader->GetOutput();
        image->Register(NULL);
        //image->SetSource(NULL);
        reader->GetOrigin(xStart, yStart);
        reader->GetStep(xStep, yStep);
        reader->Delete();
        if (xStep < 0.)
        {
            int dims[3];
            image->GetDimensions(dims);
            xStart = xStart + dims[0] * xStep;
            xStep = - xStep;
        }
        if (yStep < 0.)
        {
            int dims[3];
            image->GetDimensions(dims);
            yStart = yStart + dims[1] * yStep;
            yStep = - yStep;
        }
    }
    else
        EXCEPTION1(InvalidFilesException, fname.c_str());

    // see if the selection is such that we read the whole image anyways
    if (haveSelections && (xmin == 0) && (ymin == 0))
    {
        int extents[6];
        image->GetExtent(extents);
        if ((xmax > extents[1]) && (ymax > extents[3]))
            haveReadWholeImage = true;
    }

#endif
}


// ****************************************************************************
//  Method: avtImageFileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it.  This method is the mechanism for doing
//      that.
//
//  Programmer: Chris Wojtan
//  Creation:   Thu Jun 3 09:50:31 PDT 2004
//
//  Modifications:
//    Mark C. Miller, Tue Nov  9 13:41:33 PST 2004
//    Removed unnused pointvarnames and pointvars
//
// ****************************************************************************

void
avtImageFileFormat::FreeUpResources(void)
{
    cellvars.clear();

    if (image != NULL)
    {
        image->Delete();
        image = NULL;
    }

    haveReadWholeImage = false;
}


// ****************************************************************************
//  Method: avtImageFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: Chris Wojtan
//  Creation:   Thu Jun 3 09:50:31 PDT 2004
//
//  Modifications:
//  
//     Chris Wojtan Mon Jun 21 15:39 PDT 2004
//     Moved variable name loading from constructor to this function,
//     Changed the 3D vector to 4D vector
//
//     Mark C. Miller, Tue Nov  9 13:41:33 PST 2004
//     Moved code to populate variable names to constructor
//
//     Mark C. Miller, Wed Dec 15 10:26:07 PST 2004
//     Added support for node-centered representation of image
//
//     Hank Childs, Fri Mar 18 10:22:20 PST 2005
//     Make sure to call initialize.  Also add support for image volumes.
//
//     Jeremy Meredith, Thu Aug  7 16:12:17 EDT 2008
//     Use const char* for string literals.
//
// ****************************************************************************

void
avtImageFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md)
{
    Initialize();

    int dim = (fext == "imgvol" ? 3 : 2);
    AddMeshToMetaData(md, "ImageMesh", AVT_RECTILINEAR_MESH, NULL, 1, 0, 
                      dim, dim);
    AddMeshToMetaData(md, "ImageMesh_nodal", AVT_RECTILINEAR_MESH, NULL, 
                      1, 0, dim, dim);

    // add these variable names to the metadata
    for(size_t i=0; i<cellvarnames.size(); i++)
    {
        AddScalarVarToMetaData(md, cellvarnames[i], "ImageMesh", AVT_ZONECENT);
        AddScalarVarToMetaData(md, cellvarnames[i] + "_nodal", 
                               "ImageMesh_nodal", AVT_NODECENT);
    }

    // Check to see if we support colors.
    bool supportColors = true;
    int numExtsThatDontSupportColors = 5;
    const char *extsThatDontSupportColors[5] = { "spr", "SPR", "sdt", "SDT",
                                                 "imgvol" };
    for (int i = 0 ; i < numExtsThatDontSupportColors ; i++)
        if (image_fext == extsThatDontSupportColors[i])
            supportColors = false;

    // we do not create a list of vector variables,
    // because there is only one vector variable, "color", 
    // which is created from these scalar variables
    if (supportColors)
    {
        AddVectorVarToMetaData(md, "color", "ImageMesh", AVT_ZONECENT, 4);
        AddVectorVarToMetaData(md, "color_nodal", "ImageMesh_nodal",
                               AVT_NODECENT, 4);
    }
#ifdef EXPOSE_DEPTH
    AddScalarVarToMetaData(md, "depth", "ImageMesh", AVT_ZONECENT);
    AddScalarVarToMetaData(md, "normalized_depth", "ImageMesh", AVT_ZONECENT);
#endif
    if (fext == "imgvol")
    {
        md->SetFormatCanDoDomainDecomposition(true);
    }
}


// ****************************************************************************
//  Method: avtImageFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  
//
//  Arguments:
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: Hank Childs
//  Creation:   March 18, 2005
//
// ****************************************************************************

vtkDataSet *
avtImageFileFormat::GetMesh(const char *meshname)
{
    if (haveImageVolume)
    {
        return GetImageVolumeMesh(meshname);
    }
    else
    {
        return GetOneMesh(meshname);
    }
}


// ****************************************************************************
//  Method: avtImageFileFormat::GetImageVolumeMesh
//
//  Purpose:
//      Gets a series of meshes.
//
//  Programmer: Hank Childs
//  Creation:   March 23, 2005
//
//  Modifications:
//
//    Hank Childs, Tue Jun  7 11:45:23 PDT 2005
//    Added ghost nodes.
//
//    Mark C. Miller, Fri Apr 23 17:27:41 PDT 2010
//    Replaced ghost nodes with ghost zones.
//
//    Cyrus Harrison, Thu Oct  6 10:31:59 PDT 2011
//    Removed zonal shift of -.5, added ZStart offset.
//
// ****************************************************************************

vtkDataSet *
avtImageFileFormat::GetImageVolumeMesh(const char *meshname)
{
    int i;
    int rank = PAR_Rank();
    int nprocs = PAR_Size();
    bool doGhostZones = nprocs>1;
    bool isNodal = !strcmp(meshname, "ImageMesh_nodal");
    int globalZoneCount = (int)subImages.size();
    if (isNodal) globalZoneCount--; 
    //int globalZoneStart = 0;
    float globalZStart = (specifiedZStart ? zStart : 0.);
    float ZStep = (specifiedZStep ? zStep : 1.);

    int localZoneStart;
    int localZoneCount = globalZoneCount / nprocs;
    int procsWithExtra = globalZoneCount % nprocs;

    if (rank < procsWithExtra)
    {
        localZoneCount++; // this proc has extra layer of zones
        localZoneStart = rank * localZoneCount;
    }
    else
    {
        localZoneStart = procsWithExtra * (localZoneCount+1) +      // procs w/ extra layer
                         (rank - procsWithExtra) * localZoneCount;  // procs w/o extra layer
    }

    // -1 means we will take any image.
    indexOfImageToRead = -1;
    vtkDataSet *one_slice = GetOneMesh(meshname);

    //
    // This shouldn't happen, but it is always good to check assumptions.
    //
    if (one_slice == NULL)
    {
        debug1 << "Return value from avtImageFileFormat::GetOneMesh was "
               << "NULL, returning early" << endl;
        return NULL;
    }
    if (one_slice->GetDataObjectType() != VTK_RECTILINEAR_GRID)
    {
        debug1 << "Return value from avtImageFileFormat::GetOneMesh was "
               << "not rectilinear, returning early" << endl;
        return NULL;
    }

    vtkRectilinearGrid *rgrid = (vtkRectilinearGrid *) one_slice;
    int dims[3];
    rgrid->GetDimensions(dims);
    dims[2] = localZoneCount+1;
    int start = localZoneStart;
    if (localZoneStart > 0) {dims[2]++; start--;}
    if (localZoneStart + localZoneCount < globalZoneCount-1) dims[2]++;
    vtkFloatArray *z = vtkFloatArray::New();
    z->SetNumberOfTuples(dims[2]);
    for (i = 0 ; i < dims[2]; i++)
        z->SetTuple1(i, globalZStart + (start + i) * ZStep);
    rgrid->SetDimensions(dims);
    rgrid->SetZCoordinates(z);
    z->Delete();

    if (doGhostZones)
    {
        int nvals = (dims[0]-1)*(dims[1]-1)*(dims[2]-1);
        int nplane = (dims[0]-1)*(dims[1]-1);
        vtkUnsignedCharArray *ghost_zones = vtkUnsignedCharArray::New();
        ghost_zones->SetName("avtGhostZones");
        ghost_zones->SetNumberOfTuples(nvals);
        unsigned char *gzp = ghost_zones->GetPointer(0);
        for (i = 0 ; i < nvals ; i++)
            gzp[i] = 0;
        unsigned char val = 0;
        avtGhostData::AddGhostZoneType(val, DUPLICATED_ZONE_INTERNAL_TO_PROBLEM);
        if (localZoneStart > 0)
        {
            for (i = 0 ; i < nplane ; i++)
                gzp[i] = val;
        }
        if (localZoneStart + localZoneCount < globalZoneCount-1)
        {
            int planeStart = nvals-nplane;
            for (i = planeStart ; i < nvals ; i++)
                gzp[i] = val;
        }
        rgrid->GetCellData()->AddArray(ghost_zones);
        ghost_zones->Delete();
    }

    return rgrid;
}


// ****************************************************************************
//  Method: avtImageFileFormat::GetOneMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: Chris Wojtan
//  Creation:   Thu Jun 3 09:50:31 PDT 2004
//
//  Modifications:
//    Mark C. Miller, Thu Nov  4 18:21:25 PST 2004
//    I made it get xdim/ydim (and mins/maxs) from vtkImageData object
//
//    Mark C. Miller, Wed Dec 15 10:26:07 PST 2004
//    Added support for node-centered representation of image
//
//    Hank Childs, Fri Mar 18 11:41:04 PST 2005
//    Renamed to GetOneMesh.  Also added support for setting up images that
//    do not start at (0,0) with steps of size 1.
//
//    Mark C. Miller, Fri Apr 23 23:30:59 PDT 2010
//    Changed logic to use 'isNodal' so its consistent with GetImageVolumeMesh
//
//    Cyrus Harrison, Thu Oct  6 10:31:59 PDT 2011
//    Removed zonal shift of -.5.
//
// ****************************************************************************

vtkDataSet *
avtImageFileFormat::GetOneMesh(const char *meshname)
{
    ReadInImage();

    bool isNodal = !strcmp(meshname, "ImageMesh_nodal");
    int dims[3];
    image->GetDimensions(dims);
    int xdim = dims[0];
    int ydim = dims[1];

    // Set up rectilinear grid representing the image...
    // Since the colors are cell variables,
    //    we create an extra row of nodes in each dimesion
    //    so we can have the correct number of cells.
    int i;
    vtkFloatArray *xCoords = vtkFloatArray::New();
    for(i=0; i<xdim + !isNodal; i++)
        xCoords->InsertNextValue((float) xStart + i*xStep);
    vtkFloatArray *yCoords = vtkFloatArray::New();
    for(i=0; i<ydim + !isNodal; i++)
        yCoords->InsertNextValue((float) yStart + i*yStep);
    vtkFloatArray *zCoords = vtkFloatArray::New();
    zCoords->InsertNextValue(0.0);

    vtkRectilinearGrid *dataset = vtkRectilinearGrid::New();
    dataset->SetDimensions(xdim+!isNodal,ydim+!isNodal,1);
    dataset->SetXCoordinates(xCoords);
    dataset->SetYCoordinates(yCoords);
    dataset->SetZCoordinates(zCoords);

    xCoords->Delete();
    yCoords->Delete();
    zCoords->Delete();

    return dataset;
}


// ****************************************************************************
//  Method: avtImageFileFormat::GetVar
//
//  Purpose:
//      Gets the var associated with this file.  
//
//  Arguments:
//      varname    The name of the var of interest.  This can be ignored if
//                 there is only one var.
//
//  Programmer: Hank Childs
//  Creation:   March 18, 2005
//
// ****************************************************************************

vtkDataArray *
avtImageFileFormat::GetVar(const char *varname)
{
    if (haveImageVolume)
    {
        return GetImageVolumeVar(varname);
    }
    else
    {
        return GetOneVar(varname);
    }
}


// ****************************************************************************
//  Method: avtImageFileFormat::GetImageVolumeVar
//
//  Purpose:
//      Gets the variable for the image volume.
//
//  Programmer: Hank Childs
//  Creation:   March 23, 2005
//
//  Modifications:
//    Mark C. Miller, Fri Apr 23 17:28:15 PDT 2010
//    Added support for ghost zones.
//
//    Brad Whitlock, Thu Sep 29 15:39:56 PDT 2011
//    I added support for a 4-component color vector.
//
// ****************************************************************************

vtkDataArray *
avtImageFileFormat::GetImageVolumeVar(const char *varname)
{
    int rank = PAR_Rank();
    int nprocs = PAR_Size();

    // Number of components.
    int nc = (strncmp(varname, "color", 5) == 0) ? 4 : 1;

    // examine name to see if its node-/zone-centered
    const char *vtmp = varname;
    size_t len = strlen(varname);
    size_t len2 = strlen("_nodal");
    if (len > len2)
        vtmp = varname + len - len2;
    bool isNodal = !strcmp(vtmp, "_nodal");

    int globalZoneCount = (int)subImages.size();
    if (isNodal) globalZoneCount--; 
    //int globalZoneStart = 0;
    //float globalZStart = (specifiedZStart ? zStart : 0.);
    //float ZStep = (specifiedZStep ? zStep : 1.);

    int localZoneStart;
    int localZoneCount = globalZoneCount / nprocs;
    int procsWithExtra = globalZoneCount % nprocs;

    if (rank < procsWithExtra)
    {
        localZoneCount++; // this proc has extra layer of zones
        localZoneStart = rank * localZoneCount;
    }
    else
    {
        localZoneStart = procsWithExtra * (localZoneCount+1) +      // procs w/ extra layer
                         (rank - procsWithExtra) * localZoneCount;  // procs w/o extra layer
    }

    int count = localZoneCount+isNodal;
    int start = localZoneStart;
    if (localZoneStart > 0) {count++; start--;}
    if (localZoneStart + localZoneCount < globalZoneCount-1) count++;

    vtkFloatArray *arr = vtkFloatArray::New();
    arr->SetNumberOfComponents(nc);
    bool haveInitialized = false;
    int  valsPerSlice = 0;
    for (int i = start; i < start+count; i++)
    {
        indexOfImageToRead = i;
        vtkDataArray *one_slice = GetOneVar(varname);

        //
        // This shouldn't happen, but it is always good to check assumptions.
        //
        if (one_slice == NULL)
        {
            debug1 << "Return value from avtImageFileFormat::GetOneVar was "
                   << "NULL, returning early" << endl;
            return NULL;
        }
        if (one_slice->GetDataType() != VTK_FLOAT)
        {
            debug1 << "Return value from avtImageFileFormat::GetOneVar was "
                   << "not floating point, returning early" << endl;
            return NULL;
        }

        if (!haveInitialized)
        {
            valsPerSlice = one_slice->GetNumberOfTuples() * nc;
            arr->SetNumberOfTuples(one_slice->GetNumberOfTuples() * count);
            haveInitialized = true;
        }
       
        float *p1 = (float *) one_slice->GetVoidPointer(0);
        float *p2 = (float *) arr->GetVoidPointer(0);
        p2 += valsPerSlice*(i-start);
        memcpy(p2, p1, valsPerSlice*sizeof(float));
        one_slice->Delete();
    }

    return arr;
}


// ****************************************************************************
//  Method: avtImageFileFormat::GetOneVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      varname    The name of the variable requested.
//
//  Programmer: Chris Wojtan
//  Creation:   Thu Jun 3 09:50:31 PDT 2004
//
//  Modifications:
//    Mark C. Miller, Thu Nov  4 18:21:25 PST 2004
//    I made it get xdim/ydim from vtkImageData object
//
//    Mark C. Miller, Tue Nov  9 13:41:33 PST 2004
//    Added code to return float data directly from vtkImageData object
//    instead of through intermediary float vectors
//
//    Mark C. Miller, Wed Dec 15 10:26:07 PST 2004
//    Added support for node-centered representation of image
//
//    Hank Childs, Wed Mar 23 15:30:41 PST 2005
//    Renamed to GetOneVar from GetVar.
//
//    Hank Childs, Fri May 20 11:03:10 PDT 2005
//    Fixed flipping of green and blue channels.  Also added fast track for
//    float and unsigned char data.
//
//    Brad Whitlock, Thu Sep 29 15:39:14 PDT 2011
//    I added support for returning a 4-component color vector.
//
//    Brad Whitlock, Wed Feb  7 16:48:09 PST 2018
//    Read Z buffer if requested.
//
// ****************************************************************************

vtkDataArray *
avtImageFileFormat::GetOneVar(const char *varname)
{
    ReadInImage();

    int dims[3];
    image->GetDimensions(dims);
    int xdim = dims[0];
    int ydim = dims[1];

    int extents[6];
    image->GetExtent(extents);
    int xmin = extents[0];
    //int xmax = extents[1];
    int ymin = extents[2];
    //int ymax = extents[3];

#ifdef EXPOSE_DEPTH
    if(strcmp(varname, "depth") == 0 || strcmp(varname, "normalized_depth") == 0)
    {
        bool failed = true;
        vtkFloatArray *farr = vtkFloatArray::New();
        farr->SetNumberOfTuples(dims[0]*dims[1]);

        std::string zfile(fname);
        zfile = zfile.substr(0, zfile.size()-image_fext.size()-1) + std::string(".depth.Z");
        FileFunctions::VisItStat_t props;
        if(FileFunctions::VisItStat(zfile, &props) == 0)
        {
            vtkZLibDataCompressor *compressor = vtkZLibDataCompressor::New();
            if(compressor != NULL)
            {
                FILE *f = fopen(zfile.c_str(), "rb");
                if(f != NULL)
                {
                    unsigned char *cd = new unsigned char[props.st_size];
                    fread(cd, props.st_size, sizeof(unsigned char), f);
                    fclose(f);

                    compressor->Uncompress(cd, props.st_size,
                        (unsigned char *)farr->GetVoidPointer(0),
                        dims[0]*dims[1]*sizeof(float));

                    float *tmp = new float[dims[0]];
                    for(int j = 0; j < dims[1]/2; ++j)
                    {
                        float *a = ((float*)farr->GetVoidPointer(0)) + j*dims[0];
                        float *b = ((float*)farr->GetVoidPointer(0)) + (dims[1]-1-j)*dims[0];
                        memcpy(tmp, a, dims[0]*sizeof(float));
                        memcpy(a, b, dims[0]*sizeof(float));
                        memcpy(b, tmp, dims[0]*sizeof(float));
                    }

                    if(strcmp(varname, "normalized_depth") == 0)
                    {
                        int n = dims[0]*dims[1];
                        float *fptr = (float *)farr->GetVoidPointer(0);
                        for(int i = 0; i < n; ++i)
                        {
                            if(fptr[i] >= 256.f)
                                fptr[i] = 1.f;
                            else if(fptr[i] < 0.f)
                                fptr[i] = 0.f;
                            else
                                fptr[i] = fptr[i] / 256.f;
                        }
                    }

                    compressor->Delete();
                    delete [] cd;
                    delete [] tmp;
                    failed = false;
                }
            }
        }

        if(failed)
        {
            float *fptr = (float *)farr->GetVoidPointer(0);
            int n = dims[0]*dims[1];
            for(int i = 0; i < n; ++i)
                fptr[i] = 1.f;
        }

        return farr;
    }
#endif

    int nChannels = image->GetNumberOfScalarComponents();

    int nc = 1;
    int channel = 0;
    if (strncmp(varname, "color", 5) == 0)
    {
        nc = 4;
        channel = -2;
    }
    else 
    {
        if (strncmp(varname, "red", 3) == 0)
            channel = 0;
        else if (strncmp(varname, "green", 5) == 0)
            channel = 1;
        else if (strncmp(varname, "blue", 4) == 0)
            channel = 2;
        else if (strncmp(varname, "alpha", 5) == 0)
            channel = 3;
        else if (strncmp(varname, "intensity", 9) == 0)
        {
            // Use red as intensity if there are not 3 channels in the image.
            channel = (nChannels < 3) ? 0 : -1;
        }
        else
        {
            EXCEPTION1(InvalidVariableException, varname);
        }

        // If the image is just 1 channel then use red.
        if (nChannels == 1)
            channel = 0;
    }

    vtkFloatArray *scalars = vtkFloatArray::New();
    scalars->SetNumberOfComponents(nc);
    scalars->SetNumberOfTuples((xdim)*(ydim));
    float *ptr = (float *)scalars->GetVoidPointer(0);

    int i, j;

    if (channel >= 0)
    {
        // Extract a component.
        if (image->GetScalarType() == VTK_FLOAT)
        {
            float *data = (float *) image->GetScalarPointer();
            for (j = 0; j < ydim; j++)
            {
                for (i = 0; i < xdim; i++)
                {
                    int index = j*xdim + i;
                    ptr[index] = data[nChannels*index + channel];
                }
            }
        }
        else if (image->GetScalarType() == VTK_UNSIGNED_CHAR)
        {
            unsigned char *data = (unsigned char *) image->GetScalarPointer();
            for (j = 0; j < ydim; j++)
            {
                for (i = 0; i < xdim; i++)
                {
                    int index = j*xdim + i;
                    ptr[index] = (float) (data[nChannels*index + channel]);
                }
            }
        }
        else
        {
            for (j = 0; j < ydim; j++)
            {
                for (i = 0; i < xdim; i++)
                    ptr[j*xdim + i] = 
                     image->GetScalarComponentAsDouble(i+xmin,j+ymin,0,channel);
            }
        }
    }
    else if(channel == -1)
    {
        // Intensity.
        if (image->GetScalarType() == VTK_FLOAT)
        {
            float *data = (float *) image->GetScalarPointer();
            for (j = 0; j < ydim; j++)
            {
                for (i = 0; i < xdim; i++)
                {
                    int index = j*xdim + i;
                    float r = data[nChannels*index + 0];
                    float g = data[nChannels*index + 1];
                    float b = data[nChannels*index + 2];
                    ptr[j*xdim + i] = (r + g + b) / 3.0;
                }
            }
        }
        if (image->GetScalarType() == VTK_UNSIGNED_CHAR)
        {
            unsigned char *data = (unsigned char *) image->GetScalarPointer();
            for (j = 0; j < ydim; j++)
            {
                for (i = 0; i < xdim; i++)
                {
                    int index = j*xdim + i;
                    unsigned char r = data[nChannels*index + 0];
                    unsigned char g = data[nChannels*index + 1];
                    unsigned char b = data[nChannels*index + 2];
                    ptr[j*xdim + i] = (r + g + b) / 3.0;
                }
            }
        }
        else
        {
            for (j = 0; j < ydim; j++)
            {
                for (i = 0; i < xdim; i++)
                {
                    float r, g, b;
                    r = image->GetScalarComponentAsDouble(i+xmin,j+ymin,0,0);
                    g = image->GetScalarComponentAsDouble(i+xmin,j+ymin,0,1);
                    b = image->GetScalarComponentAsDouble(i+xmin,j+ymin,0,2);
                    ptr[j*xdim + i] = (r + g + b) / 3.0;
                }
            }
        }
    }
    else if(channel == -2)
    {
        // We want a 4-component vector.
#define CONVERT_TO_COMP4(NC, BODY) \
        if(nChannels == NC) {\
            for (j = 0; j < ydim; j++) {\
                for (i = 0; i < xdim; i++) { \
                    BODY }}}
        CONVERT_TO_COMP4(1,
           *ptr++ = image->GetScalarComponentAsFloat(i+xmin,j+ymin,0, 0);
           *ptr++ = 0.f;
           *ptr++ = 0.f;
           *ptr++ = 255.f;
        );
        CONVERT_TO_COMP4(2,
           *ptr++ = image->GetScalarComponentAsFloat(i+xmin,j+ymin,0, 0);
           *ptr++ = image->GetScalarComponentAsFloat(i+xmin,j+ymin,0, 1);
           *ptr++ = 0.f;
           *ptr++ = 255.f;
        );
        CONVERT_TO_COMP4(3,
           *ptr++ = image->GetScalarComponentAsFloat(i+xmin,j+ymin,0, 0);
           *ptr++ = image->GetScalarComponentAsFloat(i+xmin,j+ymin,0, 1);
           *ptr++ = image->GetScalarComponentAsFloat(i+xmin,j+ymin,0, 2);
           *ptr++ = 255.f;
        );
        CONVERT_TO_COMP4(4,
           *ptr++ = image->GetScalarComponentAsFloat(i+xmin,j+ymin,0, 0);
           *ptr++ = image->GetScalarComponentAsFloat(i+xmin,j+ymin,0, 1);
           *ptr++ = image->GetScalarComponentAsFloat(i+xmin,j+ymin,0, 2);
           *ptr++ = image->GetScalarComponentAsFloat(i+xmin,j+ymin,0, 3);
        );
    } 

    return scalars;

}


// ****************************************************************************
//  Method: avtImageFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      varname    The name of the variable requested.
//
//  Programmer: Chris Wojtan
//  Creation:   Thu Jun 3 09:50:31 PDT 2004
//
//  Modifications:
//      Chris Wojtan Mon Jun 21 15:41 PDT 2004
//      Changed the vector dimension from 3D to 4D
//
//      Mark C. Miller, Thu Nov  4 18:21:25 PST 2004
//      I made it get xdim/ydim from vtkImageData object
//
//      Mark C. Miller, Tue Nov  9 13:41:33 PST 2004
//      Added code to return float data directly from vtkImageData object
//      instead of through intermediary float vectors
//
//      Brad Whitlock, Thu Sep 29 16:04:46 PDT 2011
//      I made it call GetVar and I beefed up that method to deal with vectors.
//
// ****************************************************************************


vtkDataArray *
avtImageFileFormat::GetVectorVar(const char *varname)
{
    return GetVar(varname);
}
